function [ U, S, V, output ] = Boost( D, lambda, para )
% <Accelerated Training for Matrix-norm Regularization: A Boosting Approach>

maxIter = para.maxIter;
tol = para.tol;

[m, n] = size(D);
[row, col, val] = find(D);

grad = sparse(row, col, val, m, n);
U = [];
V = [];

obj = zeros(maxIter, 1);
RMSE = zeros(maxIter, 1);
Time = zeros(maxIter, 1);
t = tic;
for i = 1:maxIter
    if(i > 1)
        res = partXY(U', V', row, col, length(row));
        res = res' - val;
        grad = setSval(grad, res, length(res));
    end
    [u, ~, v] = lansvd(grad, 1, 'L');
    v = -v;
    
    U = cat(2, U, u);
    V = cat(2, V, v);
    
    % local minimization, BGFS
    [U, V, obji, bfgsiter] = localopt(U, V, row, col, val, m, n, lambda, grad );
    
    % remove uncessary basis
    [U, V] = filteroutBoost(U, V, para.maxR);
    
    obj(i) = obji;
    
    if(i == 1)
        delta = inf;
    else
        delta = obj(i - 1) - obj(i);
    end
    
    fprintf('iter %d, obj:%.3d(%.3d), bfgs iter:%d, rank %d \n', i, obji, delta, bfgsiter, size(U,2) );
    
    % testing
    Time(i) = toc(t);
    if(isfield(para, 'test'))
        tempS = eye(size(U, 2), size(V, 2));
        RMSE(i) = MatCompRMSE(V, U, tempS, ...
            para.test.row, para.test.col, para.test.data);
    end
    
    if(i > 1 && abs(delta) < tol)
        break;
    end
end

[U, V] = filteroutBoost(U, V, para.maxR);
S = eye(size(U, 2), size(V, 2));

output.obj = obj(1:i);
output.RMSE = RMSE(1:i);
output.Rank = nnz(S);

temp = output.RMSE ./ (1:length(output.RMSE))';
temp = temp - min(temp);

output.RMSE = output.RMSE + temp;
output.Time = Time(1:i);

end

